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Abstract. The dynamics of a non-Markovian open quantum system described 
by a general time-local master equation is studied. The propagation of the density 
operator is constructed in terms of two processes: (i) deterministic evolution and 
(ii) evolution of a probability density functional in the projective Hilbert space. 
The analysis provides a derivation for the jump probabilities used in the recently 
developed non-Markovian quantum jump (NMQJ) method [Piilo J et ai 2008 
Phys. Rev. Lett. 100 180402]. 



PACS numbers: 03.65.Yz, 42.50.Lc 

1. Introduction 

The theory of open quantum systems [1] provides means to study in detail a selected 
part of the total Hilbert space. This is done by making a division into a system part 
and the rest of the space, which can be seen as an environment to the system of 
interest. The position of the artificial boundary between the two parts, also referred 
to as the Heisenberg cut [3], can be selected arbitrarily according to the problem in 
hand. The interaction of the system part with its environment appears as non-coherent 
dynamics, and the evolution of the reduced system has to be described by a master 
equation. 

Non-Markovian dynamics carries a trace from the past. In the open quantum 
systems, non-Markovianity emerges from the enforced reduction of state space. Non- 
Markovian dynamics has been encountered in many fields of physics, such as quantum 
optics [3], solid state physics [4], quantum chemistry [5], quantum information 
processing [6] , in the biological context [H |8] , and even as a resource to manipulate the 
quantum-classical border [9]. Solving the non-Markovian dynamics, even numerically, 
is considerably more complicated than in the Markovian case. 

Two standard types of non-Markovian master equations are (i) the Nakazima- 
Zwanzig equation [lOl [11] , which is an integro-differential equation including a time 
convolution of the state history with a memory-kernel, and (ii) time-local expressions 
\12\ I13j. In a memory- kernel master equation, the past evolution is clearly present, 
whereas in a time-local expression there is no explicit dependence on the history. 
Consequently, it appears counterintuitive that the time-local master equations could 
produce non-Markovian dynamics at all. Time-convolutionless (TCL) projection 
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operator techniques [1] 1121 113] provide the general mathematical machinery to form 
time-local master equations. 

Recently, non-Markovian quantum jump (NMQJ) method [131 [TS] was introduced 
as a simulation algorithm for the dynamics of non-Markovian open quantum systems 
which are described by a general time-local master equation. It provides the widely 
used Monte Carlo wave function (MCWF) method [16l[T7l[T8] with an extension to the 
non-Markovian regime. In contrast to other quantum-jump-based methods available 
for non-Markovian dynamics [HI [20l [21] , the NMQJ method is strictly confined to the 
original Hilbert space of the reduced system without any auxiliary degrees of freedom. 

The NMQJ method introduces the concept of a reverse quantum jump during the 
periods where decay rates, which are associated with the quantum-jump probabilities 
in the MCWF method, reach negative values. During such periods the system 
reabsorbs information and energy from its environment as decoherence and dissipation 
are reversed. Such new kind of a stochastic process occurs with a peculiar ensemble- 
dependent probability [T?l [TS]. This paper provides a detailed line of reasoning to 
derive the jump probabilities used in the NMQJ method, and introduces ensemble- 
dependent jump operators which generate the reverse jumps. 

The paper is organized as follows. In Sec. [2] the evolution of the density operator 
in terms of a pure-state decomposition is constructed. Section [3] concentrates on 
analysing the jump-like processes and shows the connection to the quantum-jump 
probabilities used in the NMQJ method. The results are discussed in Sec. HI and 
Sec. [5] summarizes and concludes the paper. 

2. Evolution 

The most general form of a time-local master equation for the reduced density operator 
p{t) achievable by the TCL procedure is [22] 




= -^[Hit),p{t)]+Y,A,it)[c,it)p{t)Clit) - l{cl{t)C,it),p{t)}), (1) 



where the coherent dynamics is generated by the Hermitian Hamiltonian H(t) 
(including possible Lamb and Stark shifts), dissipation and decoherence are induced 
by the Lindblad (jump) operators Cfc(<), and each decay channel k is equipped with 
a decay rate A.k{t) G M. The above equation is a time-dependent generalization of 
the most general form of a Markovian master equation (where > 0) generating 
completely positive dynamical maps [23l [24]. The master equation ([Ij is local in 
time, i.e., it has no explicit dependence on the past evolution of the density operator. 
However, allowing the decay rates to enter negative values produces non-Markovian 
dynamics. 

For (time-dependent) Markovian systems, where Afc(t) > always, there exists 
a piecewise deterministic process for state vectors lip) such that the density operator 
given by a pure-state decomposition 



satisfies the master equation ([T]) [25] [26]. Above, dip = DipDip* is a singular volume 
element of the Hilbert space Ti, and P[ip;t] is a time-dependent probability density 
functional on the projective Hilbert space V{T-L). The probability density P[ip]t\ is 
concentrated on the unit sphere (HV'H = 1), and it is constant for all the states within 




k 




(2) 
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P[V';t] P[Uij;t + 6t] 

Figure 1. Two-fold evolution in the projective Hilbert space 'P('H). (a) At time 
t the probability density P[ip;t] is concentrated on state l^o)- (b) Evolution over 
a small time step St propagates the initial state |0o} (open circle) continuously in 
H to state U\4>o) (along the thick grey arrow). On the other hand, the probability 
density is modified to P[Uip;t-\-5t], such that the initial weight of state \4>o) is 
redistributed between the determistically evolved U\<l>o) and other states U\(f>i) 
and U\4>2) (thin black arrows). 



the same projective ray, i.e., states that differ only by a phase factor. The probability 
distribution shows how the density operator p{t) is constructed as a statistical mixture 
of pure states = \ilj){tp\. With the above representation of the density operator, 
the average value of any physical observable described by a self-adjoint operator A is 
given by 

(A) = tY[Ap{t)] = / d^pP[^;tMA\i;) = EmA\^)], (3) 



which is the ensemble average of a quantum mechanical average {ip\A\ip). 

Describing the density operator in terms of pure states allows us to interpret the 
evolution as a two- fold process. Let us consider an infinitesimal time interval [t, t + St). 
The density operator evolves during this period as p{t) ^ p{t + St) = p{t) + Sp{t), 
where the increment 

(4) 

is given by the master equation ([T]) in the limit of 5t 0. On the other hand, the 
pure-state decomposition ^ of the density operator, being constructed of two parts, 
gives 

6pit)^ J di^{SP[^;t]p^+P[yj;t]6p^it)). (5) 

This allows an interpretation of the evolution as (i) drift of the pure states p^ i— ^ pu^ = 
+ ^Pi'i't) and (ii) drift of probability P[V';i] ^ P[Uip;t + 6t] = P[V';<] + SP[ip;t] 
among the pure states, such that /9(t + (5t) = / dtpPlUtp'jt+Stjpu-ip- Figure [1] illustrates 
the two processes involved. The task of this paper is to derive the pair of 6P[ip; t] and 
Sp^{t) such that the identity holds with the master equation ([T]). From now on, 
the time arguments are omitted for simplicity. 

2.1. Drift of pure states 

Expressing the pure-state increment as Sp^j, — \S'ip){rlj\ + \'ip){6Tjj\, the structure of the 
master equation ([1]) suggests a specific form for the evolution \tp) >-> U{t + 6t, t)\i/j) = 
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IV') + ISijj): gathering terms acting symmetrically from left and right, the evolution is 
given by U'{t + 5t,t) — 1 ~ iStH^Q/h, where the non-Hermitian Hamiltonian 

^cff = ^-yE^'=C'lC, (6) 

fc 

is familiar from the MCWF methods jlH [171 HI] • The operation of U' is non- unitary 
and does not preserve the normalization. However, this can be compensated by 
defining a propagator U as U{t + St,t)\i)) = U'{t + St,t)\i)) /\\U'{t + 6t,t)^\\, which 
gives in the first order of 6t 

|<5^) = 5t[ + i ^ A.IICVf ) IV')- (7) 

fc 

As a result, adding a non-linear term to iJ^ff confines the state to the unit sphere 
of the Hilbert space [25]. Consequently, the pure-state density operator evolves 
during an infinetesimal time increment St by 

Sp^ ^St(^-^ [H, pv] - ^ E ^k{clCk,p^}+p^ E AfcllCfc^f ) .(8) 

k fc 

It is straightforward to verify that ti[Sp^] — and (5p^ — > as (5t — >■ 0, as requested. 



2.2. Drift of probability 

The probability density drift SP can now be solved from equation ([S]) by demanding 
the identity (j4|) with the master equation ([1} , and by substituting the state increment 
([H therein. Within the pure-state decomposition ([2]), the condition becomes 



(9) 



dij(^sp[^] + stJ2P[iA^k\\CkiA\ 

k 

= StY^^kCkpCl = St j AcpY,P[<t>\^kCkP4>Cl. 



Consequently, using a delta-functional on Hilbert space, defined by / dV"5[|V-') — 
= F[4)], where F is an arbitrary smooth functional, the condition (|9]) can 
be written in form 



Aij{5P[i;]+5tY,PblA^k\\Ckn^ 

k 

-5t ( d0^P[0]Afc||Cfc0|| 



1^) 



Cfc 



d7/'Q[i/;,t,(5t]p^ = 0. 



(10) 



Here, as well as in the rest of the paper, the phase difference between the states 
appearing in the delta-functional is arbitrary and discarded. Since (fTO]) is in the form 
of the pure-state decomposition, the distribution functional Q has to vanish identically 
for all states and moments of time t. Therefore, the probability drift during the 
time interval [i, t + St) is given by 

Cfc|0 



SP[i^]^St / d(/.;^P[0]Afe||Cfc0f5 1^) 



\\CkH 



5tY,P[yj]^k\\Ck^p\ 



(11) 
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The overall probability is conserved, since J di/jSPlip] = 0, so any increase of 
the density in one region of the state space is always compensated by a decrease 
elsewhere; the connection between the regions is given by the delta-functional. The 
phase-invariance of P is conserved, and the limit 6P[tp] Q as 5t ^ guarantees the 
continuity in time. 



3. Analysis of the jump processes and their connection to NMQJ method 

The probability drift SP connects separate parts of the Hilbcrt space and, therefore, 
corresponds to a jump-like evolution. These jumps connect a source state to a target 
state. Therefore, the jump processes can be understood as a mapping Ps Pt, 
where Pg and Pt are probability density functionals for the source (S) and the 
target (T) states. Allowing all the states to make a quantum jump means that 
■Ps[V'] = ^[V'Ij where P is the momentary probability density functional describing 
the density operator in 

Since the description coincides with the earlier definition, P[ip;t] i—> _P[V-';t] + 
SP[^p;t], one can deduce that 

SP[^] = Prm - Psm. (12) 

The source and the target state probability densities can be expressed as the marginals 
of a joint probability density P^.s, such that Px[Vj] =Pr.s[V':'^] = / d0PT,s[V'j0] and 
Psli^] = -fT,s[^;V']i which gives 

dP[i;] = J d0(PT,5 [V', <P] - Pt,s V-]) • (13) 

Furthermore, by Bayes' theorem, Py 5[0, V'] — Ps[''P]Pt,s['P\'4']j which defines the 
conditional probability density functional for jumping to the target state \<p) given 
that the source state is 

Writing the probability increment (fTTj) as 

CM 



dP[^] = /"d0<5t^(p[0]Afc||Cfc0f5[|V') 



-P[V']A..||CfcVf^[|0)- pl^]) (14) 

shows that the total probability drift is indeed given by an integral over a functional 
which is antisymmetric with respect to swapping the arguments. Evidently, it is not 
at all unambiguous to form the joint probability density Pt.s based on the knowledge 
about the difference of its marginals. Adding an arbitrary antisymmetric functional 
F[ip,(j)] = /[V',</'] — f[4'T''P]j for which / d^P[?/',(/)] = 0, to the integrand of (fH)) does 
not affect the difference 5P in any way even though it would correspond to some jump 
processes. This extra term could describe, e.g., cyclic drift \ip) 10) lip), 
which would not have any net effect on any state even though pairwise net drift 
lip) •(-> 10) and, hence, Pt,s[0, V'] would be different. In any case, some non-zero / 
has to be added in order to make the connection between the properly normalized 
probability density Pt.s and the terms appearing in the integrand of (fH|) . 

The form of (|14p shows that each decay channel k contributes to the total drift 
independently. On the other hand, each term corresponding to a single channel k 
is on its own antisymmetric. The simplest possible way to proceed is to gather up 
the positive contributions by each channel together and to select /[0, V'] such that 
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it corresponds to a complementing trivial process i~> \^). It will turn out that 
this choice gives the connection to the probabilities used in the NMQJ (and MCWF) 
method. 

With the above selection, the joint probability density can be written as 



(15) 



where g is the positive contribution of each channel k and Pt.s[0, ip] c< <^[|0) — IV")] 
is the trivial process chosen such that the density Pt,s is normalized to unity. In 
the non-Markovian case the decay rates Afc may have periods of negative values. On 
the other hand, all the other components of the products appearing in the integrand 
of (fT4| (timestep St, probability density P, norm, and the delta-functional) are by 
definition non-negative, so the sign of the decay rate alone determines the sign of each 
product term. Therefore, writing the decay rates in terms of positive and negative 
parts, Afc = A+ - A^, where A^ = (|Afe| ± Afe)/2 > [H], allows gathering up the 
positive and the negative components from ((Ti| . This gives 



p|.[0,V] 



6tA+PmCk^\\'S 10) 



\\Ckyj\\ 



-StA-pmCkM^S IV) - 



c. 



(16) 



Applying Bayes' theorem defines corresponding positive conditional jump probability 
components for each channel k by identity Pt,s[0j''/'] — Ps[fp]PT,s[4'\iA = 
Y.k PsWPr.si'l^l'^] + Ps[ip]PT,s[<l>\ip]- Therefore, using the earlier identity Ps = P, 
the result is 



m 



'StA 



' \\Ck4'\ 
l|Cfe0f4lV) 



Ck 



inn 



(17) 



and the additional trivial operation, given by the normalization rule Pt^sI'HI'iI'] 
Ek Pt.sM^^ + Pr.smi'] = 1> is then 



E 



b'P^^si^'m 



sm - 1^)] 



(18) 



With this definition, also the joint probability density is automatically normalized to 
unity: Pt,s[H,7^] = 1. 

The conditional jump probabilities pT)) are for (time-dependent) Markovian case 
(A^ > and A^ = always) exactly those of the MCWF method ESI [Ml- Non- 



Markovian systems allow A^ > (during which 



0), and then the given 



conditional jump probability coincides with the ones used in the NMQJ method 
[Mj [15]. The trivial operation probability Pt,s vanishes in the integration of 
equation (fT3|) and, therefore, corresponds to not making any quantum jumps but, 
in the first order of St, deterministic evolution instead. 



4. Discussion 



Looking at a pair of states, \^) and the decay channels allow jump-like transitions 
between the two with probabilities derived in the previous section. Moreover, looking 
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at the dynamics generated by a single decay channel k, there is a pairwise symmetry 
in the derived probabilities ([Tl 



PlsM] ^ ^^Pls[>PA]^ ■ (19) 



+A 



This means that every process ji/') induced by channel k is replaced by an 

inverse process \(f>) occurring with the exactly same joint probability, in case 

the sign of the decay rate is reversed. 

In an ensemble unravelling there is a set of N states, {IV"!)}! which evolve 
stochastically according to the probabilities given in the previous section. The 
ensemble estimates the probability density functional by P[tp] ~ (1/-^) J2i ^iW ~ 
IV'^)] = il/N)J2^N,S[\^) - \cj,j)l where N, = ^{\jIj,)\\jJj,) ^ counts the number 

of ensemble members within the same projective ray as |<^j), such that = ^• 

In [15] the set is called the effective ensemble. Each state 1-01), and equivalently 

\4>j), is provided with a set of conditional jump probabilities {P,^ g[(j)\tpi]}Li{PT,s[4'\4'i]} 
corresponding to the mutually exclusive choice between using one of the decay channels 
k or the complementing trivial no-jumps option. With the given probabilities the state 
transforms to state \(f)). 

With Afc > 0, the connection between the states given by the delta-functional 
is 10) ~ Ck\'tpi)/\\Ckipi\\ (the relative phase is arbitrary). This correspondence is 
of one-to-one type in the projective space since the condition defines the ray of |0) 
unambiguously for a given \tpi). In the MCWF method the jump between the rays 
is performed by a non-linear operation H> Ck\ipi) /\\Ckipi\\ ~ |0). Equivalently, 
the same result is achieved by opearating with an ensemble-dependent linear jump 
operator Afc[0i] = Ck\'4'i){i{Ji\/\\Ckipi\\ on the source state Itjji)- 

In the non-Markovian regime with < 0, the conditional jump probability 
becomes ensemble dependent through the ratio P[0]/P['0i], and the connection 
between the states is reversed: Itpi) ~ Cfc|0)/||Cfc0||. Compared with the positive- 
decay-rate case above, the ray of the target state is now ambiguous, since the condition 
can be fulfilled by states \4>) which are not mutually equivalent and which, therefore, 
correspond to different values of P[(j)]- Consequently, there can be many target states 
10) for which the conditional jump probability for a given single decay channel k is 
non-zero. Moreover, the mapping between the states \ipi) and |0) can not be given by 
the plain Lindblad operator Ck now, but rather by the adjoint of the earlier defined 
ensemble-dependent jump operator ^[.[0] = |0) (0|C|/||Cfc0|| acting on the source 
state Itpi). Since the jump probability is proportional to P[(j)], the operators Ak[(j)] 
indeed correspond to states which are present in the ensemble. 

In conclusion, a process defined by the operator Ak[ilj] and probability P|;[0|^] — 
s\-^\''l^]\ Ak>o S^*^ replaced by its dual counterpart characterized by the adjoint 
operator Alltp] and probability [1/110] = Pp.slV'l'/'] | Afc<o' decay rate goes 

negative. This duality is illustrated in figure [2] In the non-Markovian case both the 
conditional jump probabilities and the corresponding jump operators are ensemble 
dependent. 

When looking at the drift of probability density, there has been essentially three 
different levels of detail under consideration. The overall probability drift (5P[-0] 
corresponding to the ray of {ip) is given by (|11|) . and it is for a chosen deterministic 
propagator U exact. This level is illustrated in figure|3lfaj. However, the possibility of 
using some different kind of a propagator and thereby accompanying different form of 
a probability drift is still left open. At the second level of detail the total probability 
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(a) (b) 



Figure 2. (a.) The joint probability to jump from a source state l^/i) to a target 
state |</<) using a channel k, with = A > 0, is given by P[i/']P*[0| '!/']• Operator 
Ai;[i/)] performs the jump, (h) In the case of an opposite decay rate, Aj; = —A, 
a dual process in the opposite direction is applied causing an equally strong but 
opposite probability stream P[(^]P* [^/ilt/i] . 




Figure 3. (a) The total probability drift 5P[^] corresponding to the projective 
ray of li/i) during a time interval [t, t + 5t) is exact for a given deterministic 
propagator, (b) For a pair of states, \ip) and the net exchange consists 
of two opposing fluxes of probability Py 5 [</>,?/'] and P-p <^]. The form of 
Pt,S is ambiguous for a given 5P. (c) The joint probability Pt,s divides into 
mutually exclusive channels k corresponding to the set of Lindblad operators 
in the master equation. Each channel provides the state with a conditional 
jump probability and a jump operator. Channels with positive (solid arrows) 
and negative (dashed arrows) decay rates have different but related probability- 
operator pairs (cf. figure [2]l. 



drift (5P[^] is expressed as a difference of pairwise processes between states and 
with fluxes given by the joint probabilities Pt,5[0, tp] and Pt^sH'^ 0] [see equation (fT3|l 
and figure [Slfbj] . This division was shown to be ambiguous and to allow arbitrary 
processes with zero net effect but non-zero pairwise net flux. Neglecting the possible 
non-trivial processes with vanishing net effect and dividing the overall probability flux 
Pt,s'[^, ■0] from the source state 1-0) to the target state \(j>) into separate terms, which 
correspond to the set of decay channels in the master equation ([T]), is the third level 
of detail [see equation (fT5|) and figure [3lfcJ] . Especially, breaking the joint probability 
concerning the whole ensemble by Bayes' theorem to conditional jump probabilities 
provides each ensemble member with a set of jump probabilities. This is the level 
where the MCWF and NMQJ methods operate. 

As a last remark, let us discuss the non-trivial question concerning positivity in 
the non-Markovian dynamics. According to the proof given in |22j . the violation of 
positivity, which manifests the breakdown of the approximation scheme used during 
the derivation of the master equation, necessarily reflects to a singularity in the 
conditional reverse-jump probability. In the terminology of this paper, at the moment 
the density operator becomes negative, there are necessarily a channel fcg and states 
[■0) and 10), such that A^^ < and P^"5[0,0] > even though P5[0] = 0. This is in 
conflict with Bayes' theorem, and the algorithm terminates. However, the breakdown 
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of the NMQJ method at the negativity is only an imphcation and the reverse does 
not need to be true. 

The key element is the assumption that the algorithm always produces a proper 
probability density functional P within the set of basis states generated by the 
algorithm itself (the support of the probability density functional, {IV") G 7^ | /'[■(/'] ^ 
0}, is understood as the basis set since the density operator is expressed as a 
statistical mixture of such states). Given a pure initial state p(0) = |V'o)(V'o|, 
corresponding to a point spectrum P[V';0] = S[\ip) — Hjq)], the basis set at later 
times consists of all the states reachable via different combinations of operations 
by propagator U and jump operators Ck at any consecutive times on state \^po)^■ 
U{t,tn)Ck„{tn)U{tn,tn^i) ■ ■ ■ Cki{ti)U{ti,0)\'iljo). Naturally, these states do not 
generally correspond to the momentary eigenvectors of the density operator, and hence 
the weights of such states are not necessarily non-negative even though p would be 
positive. Therefore, the question arises, whether it could be possible to generate the 
dynamics in terms of quasiprobability densities, for which P may have also negative 
values. Clearly, equations ([8]) and (fTTj) give the increments that could be used for 
direct integration within the projective Hilbert space without the estimation by a 
stochastic process. If the dynamics is confined to a restricted number of different pure 
states (cf. examples in [TS]) this method would be also efficient, while in the general 
case a stochastic process would be needed to restrict the number of different states. 

5. Summary and conclusions 

The dynamics of an open quantum system described by a general time-local non- 
Markovian master equation was analyzed. Using a pure-state decomposition the 
dynamics was described as (i) deterministic evolution and (ii) evolution of the 
corresponding probability density functional in the projective Hilbert space. This 
division allows the interpretation of the dynamics as a continuous drift in Hilbert 
space accompanied with jump-like transitions between different parts of the state 
space. 

Certain assumptions and choices were made in order to recover the NMQJ 
method: (i) the deterministic evolution is given by a certain specific propagator, 
(ii) probability density functional describing the density operator maps to a valid 
probability density functional during every time step, (iii) each decay channel k 
contributes independently and the choice between the channels is mutually exclusive, 
and (iv) the remaining jump probability is complemented by a trivial identity 
operation corresponding to a no-jump option. The first choice is motivated by the 
symmetry issues and it coincides with the one used in the MCWF method, the 
extension of which is the NMQJ method. For the used deterministic propagator, the 
corresponding probability drift is exact. Using the second assumption, the probability 
drift is given by a difference between two marginals of a joint probability distribution 
for source and target states. The third and the fourth assumptions are motivated 
purely by simplicity. 

The reverse jumps |141 I15j were formulated in terms of ensemble-dependent 
jump operators. In contrast to the Markovian case, both the conditional jump 
probabilities and the corresponding jump operators are in the non-Markovian case 
ensemble dependent, which can be seen as a manifestation of the memory effects. 

Open questions remain for future investigations. The NMQJ method terminates 
as soon as the conditional jump probability diverges irrespective of whether the density 
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operator is actually about to lose its positivity or not. Could some of the listed 
assumptions perhaps be relaxed in order to allow a more general stochastic simulation? 
Or more generally, is there any stochastic process within the reduced system's Hilbert 
space that is able to simulate the general time-local master equation whenever it is 
physically valid? 
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